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Abstract 

This arXiv upload is to clarify that the now well-known sorted QR MIMO 
decoder was first presented in the 1995 IUGG General Assembly. We clearly 
go much further in the sense that we directly incorporated reduction into this 
one step, non-exact suboptimal integer solution. Except for these first few lines 
up to this point, this paper is an unaltered version of the paper presented at 
the IUGG1995 Assembly in Boulder. 

The Ambiguity resolution of GPS carrier phase observables is crucial in high 
precision geodetic positioning and navigation applications. It consists of two 
aspects: estimating the integer ambiguities in the mixed integer observation 
model and examining whether they are sufficiently accurate to be fixed as 
known nonrandom integers. We shall discuss the first point in this paper from 
the point of view of integer programming. A one-step nonexact approach 
is proposed by employing minimum diagonal pivoting Gaussian decomposi- 
tions, which may be thought of as an improvement of the simple rounding-off 
method, since the weights and correlations of the Boating-estimated ambi- 
guities are fully taken into account. The second approach is to reformulate 
the mixed integer least squares problem into the standard 0-1 linear integer 
programming model, which can then be solved by using, for instance, the 
practically robust and efficient simplex algorithm for linear integer program- 
ming. It is exact, if proper bounds for the ambiguities are given. Theoretical 
results on decorrelation by unimodular transformation are given in the form 
of a theorem. 

1 Introduction 

Three types of observables may be derived from tracking GPS satellites: pseudorange (code) mea- 
surements, raw Doppler shifts (or equivalently range rates) and carrier phases. They are used 
at different levels of accuracy for different purposes of applications (see e.g Wells et al. 1986; 
Leick 1990; Hofmann-Wellenhof et al. 1992; Seeber 1993; Melbourne 1985). The carrier phase 
measurements, together with the accurate code observables (if available), have been dominating 
in high precision geodetic positioning and navigation applications. The mathematical model can 
symbolically be written below 



Here R and are respectively the observables of pseudoranges and carrier phases, £r and e$ 
are the random errors of the observables, X is the coordinate vector to be estimated, and ffi(-) 
and f$(.) are nonlinear functionals of X. B/j, B$ and are the coefficient matrices. A is the 
vector of nuisance parameters such as the synchronization errors of receiver and satellite clocks and 
ionospheric corrections. If overparametrization occurs to A, it is generally not estimable (Wells et 
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al. 1987). Thus we shall assume that proper reparametrization has been made by, for instance, 
choosing proper datum parameters (Wells et al. 1987) or using differencing and nuisance parameter 
elimination techniques (see e.g. Goad 1985; Schaffrin & Grafarend 1986), to ensure that the 
remaining nuisance parameters are estimable. Z is the vector of integral ambiguities inhered in the 
carrier phase observables. 

Accurate and reliable resolution of the integral ambiguity vector has been playing a crucial role 
in high precision positioning. There are currently many approximate proposals available to resolve 
Z. They may be treated in two categories: simple (sequential) rounding-off of a real number 
to its nearest integer with and/or without using constraint criteria (Blewitt 1989; Talbot 1991; 
Hwang 1991; Seeber 1993; Hofmann-Wellenhof et al. 1992), and searching methods by employing 
the information on the prior statistics and geometry (nonlinear functionals and design matrices) 
of the observables (Counselman et al. 1981; Remondi 1990, 1991; Frei &: Beutler 1990; Mader 
1990; Mervart et al. 1994). Betti, Crespi &: Sanso (1993) recently proposed a Bayesian approach 
to resolution of ambiguity. Chen & Lachapelle (1994) proposed a fast ambiguity search filtering 
approach to reducing the number of possible candidates in the searching area. It may be worth 
noting that the fast rapid ambiguity resolution method proposed by Frei &: Beutler seems to have 
enjoyed its wide approval. A key element of the method is the use of some formal statistics to pick 
up a solution. It may be proved that the statistic used for selecting the candidates of ambiguities is 
not mathematically rigorous, since the ambiguity-free and ambiguity-fixed solution vectors are both 
derived by using the same set of carrier phase observations. The method seems quite successful in 
practice, however. 

Recent progress in resolving the integral ambiguity vector has been made by Teunissen (1994). 
His approach consists of three steps: (1) decorrelation of the floating-estimated ambiguities by 
Gaussian transformation, which may be said to characterize the novelty of the new approach, (2) 
searching for the solution to the transformed integer least squares problem within a superellipsoid 
corresponding to a certain level of confidence, and (3) back-substituting the solution just derived 
for the ambiguity vector in the original model. The success of the approach will depend, to a great 
extent, on the first two steps. Testing results of the approach can be found in Teunissen (1994) 
and de Jonge & Tiberius (1994). Decorrelation techniques may be also well suited to explain an 
important finding by Melbourne (1985), that the widelane ambiguity is easier to solve, based on 
the one epoch dual frequency carrier phase and code-derived pseudorange model. 

The purpose of this paper is to further study the GPS ambiguity resolution as a mixed integer 
least squares (LS) mathematical programming problem. Unimodular integer transformation is 
used to statistically decorrelate the floating-estimated ambiguities, which summarizes the first 
two conditions of transformation proposed by Teunissen (1994). Two methods for solving the 
transformed integer LS problem are then proposed. The first one is to decompose the transformed 
positive definite matrix into a lower and an upper triangle by choosing the minimum diagonal 
elements. In this way, we are sure that a wrongly selected ambiguity will be penalized. No iterations 
are required, thus it should improve the sum of square of the residuals derived by rounding the 
transformed real values to their nearest integers. The second one is to reformulate the transformed 
integer LS problem to a quadratic 0-1 nonlinear programming, and then further to a 0-1 linear 
integer programming. Thus simplex algorithms can be employed to efficiently solve the linear 
integer programming problem, with which one need not test every point in the feasible solution set. 

2 Integer and mixed integer least squares models 

In the application of the GPS system to high precision positioning and navigation, the GPS satellites 
have been treated as space targets with known positions, unless the determination of the satellite 
orbits is of interest. In this paper, we assume that the coordinates of the satellites are given, 
which can be computed, for instance, from the (precision) ephemerides. Furthermore, given a set 
of approximate coordinates of the stations, we can linearize the observation equations (la) and (lb) 
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as 



where 



y R = A R AX + B R AX + e R 
y$ = A$AX + B$AA + B Z AZ + e$ 

y R = R - fR(Xo) - B^Ao 
y$ = * — f$(X ) — B$Ao — B^Z 
AX = X-X ; AA = A - A 
AZ = Z - Z . 



(2a) 
(26) 

(3a) 
(36) 
(3c) 
(3d) 



Xo and Ao are the approximate values of X and A, respectively. Zo are integer approximate values 
of Z, and thus AZ remain integral. 

Rewriting the linearized observation equations (2a) and (2b) in matrix form, together with the 
statistical information on the observables, we have 
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(46) 



Here Pr and P$ are respectively the weight matrices of the observables y# and y$, a 2 is the scalar 
variance component. 

Since the main interest of this paper is to discuss the mixed integer LS problem, we do not need 
to discriminate between the position unknowns X and the nuisance parameters A. Without loss 
of generality, therefore, we can simplify the model (4) as the following standard mixed real-integer 
(or simply integer in the rest of the paper) observation equations, 



y = A/3 + Bz + e 
D(y) 



P-V 



(5a) 
(56) 
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The matrices A and B are full of column rank, respectively. 
Applying the least squares criterion to (5), we have 



mm 



F = (y - A/3 - Bz) T P(y - A/3 - Bz), 



(6) 



which is the mixed integer LS problem. (6) was also called the constrained LS problem by Teunissen 
(1994). Since the variables z are discrete, we cannot use the conventional method by differentiating 
the objective function F with respect to the variables /3 and z in order to form the normal equation 
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and then solve for them. Instead, however, we differentiate F with respect to j3 and let it equal 
zero, leading to 

— = -2A T P(y-A/3-Bz) = 

or 

A T PA/3 = A T P(y - Bz). 

Hence 

/3 = (A T PA)- 1 (y-Bz). (7) 
Substituting (7) into (5) and rearranging it yield 

yi = QPBz + £i (8a) 

D[yi] = [P- 1 - A(A T PA)- 1 A T ]a 2 = Q<7 2 (86) 

where 

yi = [I - A(A T PA)- 1 A T P]y = QPy 
Q = P 1 - A(A T PA) _1 A T . 
Applying the LS method to (8), we have 

min: F x = ( yi - QPBz) T Q (yi - QPBz) 
= (y-Bz) T PQQ QP(y-Bz) 
= (y - Bz) T PQP(y - Bz) 

= y T PQPy - 2y T PQPBz + z T B T PQPBz. (9) 

The objective function F\ can further be rewritten as 

min: Fi = (z — z) T H(z — z) + y T PQ[Q - PBH _1 B T P]QPy (10) 

where 

z = H 1 B T PQPy 
H = (B T PQPB). 

Here z can readily be proved to be the floating LS estimate of the ambiguity vector AZ with 
covariance matrix H _1 cr 2 . Since y T PQ[Q — PBH x B T P]QPy is constant, the objective function 
(10) is equivalent to (Teunissen 1994; de Jonge Sz Tiberius 1994) 

min: F 2 = (z - z) T H(z - z), (11) 

which is the standard integer LS problem. 

It is now clear that the solution to the original mixed integer LS problem (6) depends solely on 
that of the standard integer LS problem (11). Denote the integer solution of z to (11) by z IN . 
Substituting it into (7), we can then obtain the LS estimates of the real parameters (3 without 
much effort. 
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3 Unimodular transformation 

In resolution of GPS carrier phase ambiguities, one of the most difficult points is to handle strong 
correlation of the matrix H. Searching for an acceptable (and/or hopefully optimal) solution of z 
is arduous, if it is solely based on the strong correlation matrix H, since testing a large number of 
combinations would have to be done. Roughly speaking, the total number of combinations required 
is computed by JJrii, where n{ is the number of integer points on an interval of line for the ith 
ambiguity, centred at the and corresponding to a significance level (see e.g Frei & Beutler 1990). 
However, if the matrix H is diagonal, one can simply round the floating values z off to the nearest 
integers, which are the integer solution of z. Therefore, an idea would emerge naturally, that one 
works with a decorrelated weight matrix instead of H. 

Such a technique was proposed recently by Teunissen (1994) (see also de Jonge & Tiberius 1994). 
His basic idea is to transform the "observables" z by G into the new ones zi = (G T z), and then 
work with the LS integer problem 

min : F 3 = (zi — zi) T Hi(zi - zi). (12) 

Here H = GH\G T . The transformation matrix G has to satisfy the following three conditions: 
(1) integer elements; (2) volume preservation; and (3) decorrelation of H into Hi. More details 
can be found in Teunissen (1994). 

Before proceeding, we shall define the unimodular matrix (see e.g. Nemhauser & Wolsey 1988). 

Definition 1. A square matrix G is said to be unimodular if it is integral and if the absolute 
value of its determinant is equal to unity, i.e. \det(G)\ = 1. 

The inverse of a unimodular matrix is also unimodular, since \det(G~ 1 )\ = l/\det(G)\ = 1, and 
because 

G 1 = G/<fei(G) = ±G. 

Here G is the adjoint matrix of G, whose elements are derived only by using the operations of 
integer multiplication, substraction and addition, and thus integer. The sign before G depends 
on the determinant of G. The second property of unimodular matrices is that the product of 
two unimodular matrices is unimodular. It is also clear that any unimodular transformation of an 
integer vector is an integer vector, too. 

By employing the concept of the unimodular matrix, we can summarize the first two conditions 
suggested by Teunissen (1994) by stating that the transformation G is unimodular. It should be 
noted that there was a misunderstanding of Teunissen's second condition of volume preservation. 
Volume preservation does not imply the preservation of the number of grid points. A simple 
example is that a unit circle centred at the origin has five grid points, while an ellipse of the same 
center with major axis 1.5 and minor axis 2/3 encloses only three grid points. 

Integer Gaussian decomposition was employed by Teunissen (1994), that indeed decorrelates the 
matrix H. What now seems to be done is to mathematically prove that we can always decorrelate 
the matrix H by using a finite number of unimodular transformations to the extent that the 
correlation coefficient of any two random variables is always less than or equal to 1/2. In order to 
do so, we need the following lemma on the inequality of matrix determinant. 

Lemma 1: For any positive definite matrix A, the following inequality 

det(A)<H aii (13) 
holds true. Here an are the diagonal elements of A. 

Proof. A positive definite matrix A can be written by Choleski decomposition as 



A = LL T 
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i-1 



where la = (an - ifi) 1 ^ 2 > 0- Thus we have 

3=1 



det(A) 



i=i 



< [[an, 



i-i 



since £ /?, > 0. □ 

3=1 

Theorem 1: For any positive definite matrix A, there exists a unimodular matrix G such that 

A = GHG T . (14) 
Here H is positive definite, too, and satisfies 



\hij\ < - min(hu, hjj) Vi,j& i^j. 



(15) 



Proof. Suppose, without loss of generality, that for any three elements an, ajj and of the 
positive definite matrix A, we have \aij\/min(aa, ajj) > 1/2. Then construct the unimodular 
matrix 

1 



Gi 



■ [dij j an]i n • • • 1 



(16a) 



if an < ajj, or 



1 • • • i a ij/ajj]i 
1 



(166) 



if ajj < an. Here [ x ]i n is the operation to round the floating number x to its nearest integer. 

Upon left- and right-multiplying A by the unimodular matrix Gi and its transpose respectively, 
the larger diagonal element is then reduced to 



vnax{an,ajj) 1\a{j j a m ni\inaij ~\~ a m i n \a,ij / a. 



2 

mmlin 



(17) 



min(aa,ajj). Repeating the same procedure to any pair of diagonal elements, we 



where a min 
have 

A n = G n ...GiAGf...G£. (18) 

Now suppose that we cannot reach the equation (14) and the inequality (15) by employing a 
finite number of unimodular matrices of the form (16), then we keep applying the same procedure 
to A n . By expression (17), it is clearly true that the minimum diagonal element of the reduced 



Presented at the IUGG95 Assembly, Boulder, Colorado, July 2-14, 1995. 



7 



matrix, say A m now, has no lower bound. It means that the minimum element can be arbitrarily 
small, which further implies by Lemma 1 that 

det(A m ) < Y[ a% < const, (19) 

where a™ are the diagonal elements of A m , const is any positive constant. Since unimodular 
transformation does preserve the determinant, we have det(A m ) = det(A) — a finite constant, 
which clearly contradicts (19). Therefore, we must be able to reach the condition (15). On the 
other hand, all the transformation matrices involved are unimodular, their product is unimodular, 
too. Denoting the final reduced matrix by H, which satisfies the condition (15), and the product 
of all the unimodular matrices by Gt, we have 

H = G t AGj (20) 

or 

A = GHG T . (21) 
Here G(= G^ -1 ) is unimodular. The proof that the matrix H is positive definite is trivial. □ 

4 Two approaches to the integer LS problem 

The integer LS problem is simply an integer quadratic programming issue. One can use any 
advanced integer programming algorithm (Parker & Rardin 1988) to solve this problem. Essentially, 
no bounds for the integer unknowns are required and no statistical techniques needed to reduce the 
number of possible candidates. More on these aspects and proper validation criteria for fixing the 
carrier phase ambiguities will be presented in a future paper. 

Though the techniques to be presented below require no decorrelation as an assumption, and 
consider that the original and the transformed LS integer problems are of the same form, the follow- 
ing discussion will be based on the transformed model, without loss of generality. After the weight 
matrix of the floating-estimated ambiguity vector is decorrelated, one can either simply round the 
transformed floating numbers off to their nearest integers, or employ searching techniques to find 
the "optimal" solution within a superellipsoid under a certain level of confidence (Teunissen 1994). 
In what follows, we shall develop two approaches to resolve the ambiguities of the transformed 
integer LS problem. 

4.1 A one-step nonexact approach by minimum diagonal pivoting Gaussian 
decomposition 

Instead of directly applying the simple rounding-off method to (12), which ignores any correlation 
information on the floating-estimated ambiguities, we propose an alternative one-step approach, 
based on the weights and correlations of the transformed ambiguities. The basic idea is to resolve 
the integer ambiguities according to their weights and correlations. As long as some of ambiguities 
are resolved, their correlations with other unfixed floating ambiguities are employed and the next 
ambiguity corresponding to the large weight is to be determined. 

In order to realize the above procedure, we have to decompose the positive definite matrix Hi 
carefully. Here we employ Gaussian decomposition by selecting the minimum diagonal element. 
The decomposition procedure consists of the following steps: 

• Selecting the minimum element among all the undecomposed diagonal elements; 

• Exchanging the rows and the columns; 

• Performing Gaussian decomposition; 
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• Replacing the square root of the decomposed element h'^,^ at the corresponding position 
of the factor matrix L; If the decomposition is not completed, then go to the first step. 
Otherwise, the decomposition is finished. 

In mathematical language, we can express the matrix Hi as 

Hi = P^LL T P£ (22) 

where is the permutation matrix which represents the exchange of the rows and columns dur- 
ing the decomposition. A significant characteristic of this decomposition is to keep the diagonal 
elements of the lower triangular matrix L in the increasing order as far as possible. 
Inserting Hi in (22) into (12), we have the objective function 

min: F 3 = (z x - zi) T P h LL T P£(zi - *i) 

= (z 2 -z 2 ) T LL T (z 2 -z 2 ) (23) 

where 

z 2 = P£zi; z 2 = P£zi. (24) 
Since the factor matrix L is lower triangular, we can rewrite (23) as 

t x t x 

min : F 4 = E E l d z 2(j) ~ z 2(j)) ? ■ (25) 

i=l j=i 

Here t z is the dimension of the ambiguity vector z (or z 2 ). The solution to the objective function 
F2 can now be derived by minimizing 

t z 

\J2h( z m -%))!> Vi - ( 26 ) 
j=i 

Hence the one-step nonexact integer ambiguity solution is immediate 



Z 2(i) 



(27) 



for all i. 

By back substituting the integer solution z^ = (z^y z 2(2)i •■•> z 2(t z )) T ' we nave the final solution 
of the integer ambiguities z, which is denoted by z IN , 

z IN = G T P h z J 2 N . (28) 



4.2 0-1 quadratic integer programming 

An obvious aim of applying the decorrelation technique to the original integer LS problem is the 
alleviation of the computational burden for finding the optimal ambiguity solution. When it is 
translated into the case of searching techniques, we expect that the total number of candidate grid 
points to be tested should be significantly reduced. Suppose that for the transformed integer LS 
problem (12) (Hi satisfies the conditions of Theorem 1), we have to search for the optimal integer 
ambiguity resolution within the hard bounds 

m° < z 1(i) < mj, V i (29) 

or in another form, 

zi(i) € \mu(= m°), m 2 i, m lsi (= m})]. (30) 



Presented at the IUGG95 Assembly, Boulder, Colorado, July 2-14, 1995. 



9 



Here Z\^ is the ith integer component of the integer vector Zi, mij, m,2i, .., and m\ Si are the 
contiguous integers — the candidate points of with the lower integer bound and the upper 
integer bound mj. Thus our mixed integer LS problem has been reduced to a quadratic integer 
programming problem with simple integer constraints. 

In what follows we shall further reformulate it by a 0-1 quadratic integer programming model. 
It has been shown by Parker & Rardin (1988) that the integer variable can be represented 
with ri 0-1 variables, i.e. 

r—l 

z 1(!) = m»+^2H l(3)l V i (31) 

3=0 

where b^ are 0-1 integer (binary) variables, r« = [log2(mj — m®)] s + 1, and [ . ] s stands for the 
integer not larger than the positive number in brackets. 

Rewriting all the integer variables in matrix form, we have 

zi=m° + Aib (32) 

where the matrix Ai is integral with elements 2 k , 

m° = (m?, ml ml) T 

b = (bf , , bJ z ) T 

Furtheron, inserting (32) into the objective function (12) yields 

min: F 3 = (A : b + m° - zi) T Hi(Aib + m° - z x ) (33) 

subject to bk = or 1 for all k. 

The objective function (33) is equivalent to 

min: F 3 = (m° - zi) T Hi(m° - z x ) + 2(m° - zi) T Hi A : b 

+b T AfHiAib. (34) 

4.3 0-1 linear integer programming 

In this subsection, we shall further reformulate the 0-1 quadratic programming (34) into a 0-1 
linear integer programming problem by using the linearization technique. The basic idea of the 
linearization technique is to introduce a new variable to replace the nonzero quadratic term bibj. 
Thus the 0-1 quadratic programming problem becomes linear. Since the new variables are obviously 
binary, all the variables in the linear programming model to be reformulated below are binary, too. 
Denoting 

v k = bibj, k = (i- l)i/2 +j, i> j 
and taking the following relations 

b1 = h 

into account, we have 

ty 

min: F 4 = (m° — zi) T Hi(m - zi) + avi (35a) 

i=l 



subject to the following constraints, 



Vi = V 1 
Vk > Vki + v k j - 1 



(356) 
(35c) 
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v k < Vki (35d) 
v k < v k j (35e) 
ki = i(i + l)/2; kj = j(j + T)/2. 
Here t v is the dimension of the 0-1 integer vector 

v = {vi,v 2 , v tv f . 
Since the first term in the objective function (35a) is constant, it is equivalent to 

min : F4 = ^ C{Vi (36a) 
i=i 

subject to the constraints (35b ~ e). (36) is obviously of the standard form of the 0-1 linear integer 
programming. It can be solved by using any standard algorithms for 0-1 linear programming 
(Pardalos & Li 1993; Nemhauser & Wolsey 1988; Parker Sz Rardin 1988; The People University of 
China 1987). However, the algorithm aspects for the program (36) will not be discussed here. 

5 Concluding remarks 

GPS carrier phase and pseudorange observables are essentially a nonlinear mixed integer observation 
model. If the GPS satellites are treated as space known targets, the model is regular. Given a 
set of approximate values of the unknown parameters such as the position coordinates and integer 
ambiguities, the nonlinear model is linearized. Estimating the parameters in the linearized mixed 
integer model is equivalent to solving a mixed integer LS problem (if the LS principle is employed), 
which can be further reduced into a standard integer LS programming. 

It has been recognized that one of the difficulties in correctly estimating the integer ambiguities 
is due to the correlations of the floating-estimated ambiguities. A decorrelation technique has 
been proposed by Teunissen (1994), based on Gaussian decomposition. We have further proved 
mathematically that there exists a unimodular matrix such that (14) and (15) hold true, which 
may be thought of as a theoretical summary (and extension) of some of the results in Teunissen 
(1994). 

Two approaches are then proposed to solve the standard linear integer LS problem (12) from the 
point of view of integer programming theory. The first approach is to Gauss-decompose the matrix 
Hi by selecting the minimum diagonal elements. In other words, we are estimating the integer 
ambiguities according to the magnitudes of the weights of the floating-estimated ambiguities and 
their correlations (as far as possible). It may be thought to be an improvement of the simple 
rounding-off method. No iterations are required. It should be noted, however, that this method 
is one-step nonexact. The extent of approximation should be further investigated. The second 
approach is to reformulate the mixed integer LS problem into a 0-1 linear integer programming 
model. Thus any standard algorithms for linear integer programming problems can be employed. 
The method will result in the exact integer solution of the ambiguities to the original mixed integer 
problem, if proper bounds for the integer unknowns in the transformed model (12) are given. 
Testing of the techniques with real data is under way. 
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